Microscopic Theory of Superconductor-Constriction-Superconductor Josephson 

Junctions in a Magnetic Field 
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Self-consistent solutions of microscopic Eilenberger theory are presented for a two-dimensional 
model of a superconducting channel with a geometric constriction. Magnetic fields, external ones 
as well as those caused by the supercurrents, are included and the relevant equations are solved nu- 
merically without further assumptions. Results concerning the influence of temperature, geometric 
parameters, of k — \l/£,o and of external magnetic fields on the Andreev bound states in the weak 
link and on the current-phase relation are presented. We find that the Andreev bound states within 
the junction obtain peculiar substructure when a finite supercurrent flows. As long as the London 
penetration depth is comparable to or bigger than the extension of the constriction, the Josephson 
effect is independent of k. Furthermore, the weak link is very insensitive to external magnetic fields. 
Features restricted to a self-consistent calculation are discussed. 

PACS numbers: 74.25.Jb, 74.45. +c, 74.50.+r 



I. INTRODUCTION 

Since their discovery, Josephson effects play a central 
role in the investigation of the superconducting state be- 
cause of their fundamental as well as technical relevance. 
Josephson effects occur whenever two extended super- 
conducting banks are weakly coupled and can be ob- 
served, among many other systems, across a geometric 
constriction of a superconductor with lateral dimensions 
much smaller than the coherence length. This kind of 
Josephson junction, usually referred to as ScS junction, 
has promising properties for highly sensitive magnetome- 
ters pj. Recent progress in microfabrication technique 
allows to pattern high quality niobium thin films with 
structures small enough for constriction-type Josephson 
junctions At the same time, a ScS Josephson junc- 

tion is a useful basis for theoretical considerations be- 
cause only the properties of the superconductor and its 
geometry are relevant. 

Superconducting point contacts serve as a simplified 
model for ScS junctions and have been investigated by 
many authors. An overview of many useful results for 
the ac and dc Josephson effect can be found in the re- 
views by Likharev [J] and Golubov et al. ||. The most 
general result was derived by Zaitsev Q and includes a 
large number of conduction channels as well as asymme- 
try and non-BCS structure of the electrodes. 
Microscopic theory for the dc Josephson effect in a model 
mesoscopic constriction has been presented in 0, HJ with 
the restriction to a small number of conducting channels 
in the constriction. This model clarified the relevance of 
Andreev bound states for the supercurrents in the junc- 
tion. Analytic solutions of microscopic theory for the dc 
Josephson effect in a ballistic superconducting microcon- 
striction have been presented in (Io| . In this approach, 
the influence of the constriction length on the critical 
current has been calculated self-consistently for temper- 
atures close to T c . 



As already stated by Kulik [ll| and Ishii [lj, 1 1 31 ] . 
a more detailed physical picture of Cooper pair tun- 
neling in Josephson junctions involves Andreev reflec- 
tion [3, EH at the interfaces between the weak link and 
the electrodes. Due to this picture, Andreev bound states 
in Josephson junctions represent electron and hole waves 
that bounce back and forth between the two electrodes 
constructively, transferring one Cooper pair in every cy- 
cle. Thus, Cooper pair transfer across Josephson junc- 
tions is mediated by Andreev bound states. With this 
view in mind, detailed knowledge of the local density of 
states (LDOS) in the vicinity of the weak link is essen- 
tial. This fact has led to the consideration of the LDOS 
in different types of Josephson junctions: am ong others, 
see 0, 1 for ScS, 0,E|, El for SNS, [H, El El E3 
for SNS and SIS and for SNS and SFS with s- and 
d-wave pairing symmetry (N = normal conductor, I = 
isolator, F = fcrromagnet). 

In the present paper, we examine a two-dimensional 
model of a ScS Josephson junction in order to investi- 
gate the dependence of the dc Josephson effect on tem- 
perature, on geometric parameters as well as the influ- 
ence of magnetic fields, external ones as well as those 
generated by the supercurrents flowing in the junction. 
Therefore, we will numerically solve microscopic Eilen- 
berger theory self-consistently for a conventional (s-wave) 
superconducting junction in the clean limit. In order to 
achieve electromagnetic self-consistency as well, we em- 
ploy a convenient method to calculate the magnetic field 
evoked by the supercurrents. This will allow us to derive 
transport properties for a realistic geometry at arbitrary 
temperature and arbitrary external magnetic field. Fur- 
thermore, the method used to solve microscopic Eilen- 
berger theory self-consistently including magnetic fields 
presented in this paper can be applied to more complex 
geometries and thus might be of more general interest. 

The paper is organized as follows. In the next sec- 
tion, we introduce our geometric model. In section IIII1 
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Figure 1: Two-dimensional model for a ScS Josephson junc- 
tion consisting of a long channel (width W = l£o throughout 
the paper) with a constriction. The constriction is defined by 
the parameters I and w. Our calculations have been carried 
out for a section of the channel with length L — 2£o (dark 
grey area) and fixed phase difference of the order parameter 
= <j)R — (j>L- The magnetic field is oriented perpendicular 
to the image plane (along the z axis). 



III. QUASICLASSICAL EILENBERGER 
THEORY 



A. Self-Consistency Equations 

Our calculations are based on quasiclassical Eilen- 
berger theory [22j which we solve self-consistently 
including the magnetic vector potential. To obtain 
Josephson weak link behaviour, it is essential to 
calculate the order parameter self-consistcntly. A 
self-consistent solution contains local modulations of 
the order parameter amplitude and can thus yield local 
suppression of superconductivity. Local suppression of 
superconductivity in turn constitutes a Josephson weak 
link. 

To solve quasiclassical Eilenberger theory self- 
consistently, we have to find solutions for the order 
parameter A(r) and the supercurrent j s (r) which satisfy 
the two self-consistency conditions: 



we present the theoretical background for our calcula- 
tions and explain how we solve the relevant equations. 
In section IIV1 we compare our results to previous ones in 
order to test our methods. In section [Vj we present rep- 
resentative results and explain how the Josephson effect 
emerges in ScS junctions. In the following, we consider 
the influence of the temperature and geometric param- 
eters (section |VI[) . of n = \l/£,o (|VH[) and of external 
magnetic fields (IVIII|) . Finally, in section IIX[ we close 
with some concluding remarks. 



A(f) = 2nVN(0)k B T £ (/(r, k F , ie n )) ps (I) 

0<£„<W C 

j s {f) = 4neN(0)k B T ^ (v F • g(r, k F , ie n )) ps (2) 

0<e n <ui<. 

These two equations, the gap equation and the current 
equation, contain the components of the quasiclassical 
propagator f{r{ x), hp, ie n ) and g(r{x), k F , ie n ) which, in 
turn, can be expressed in terms of two compl ex q uantities 
a(x) and b{x) called the Riccati amplitudes [23l |24||: 



II. MODEL GEOMETRY 

For the calculations presented in this paper, we use 
a two-dimensional model of a ScS Josephson junction 
as depicted in Fig.Q] The geometry consists of a long 
superconducting channel which is narrowed down by a 
constriction. We will calculate all relevant quantities on 
a section of the channel enclosing the constriction with 
length L = 2£ (£o = ^fF/A 00 (T = 0) is the coherence 
length). On the left and right end of this section, the 
phase of the order parameter will be set to 4>l,r- The 
width of the channel will be set to W = l£o throughout 
the paper whilst the constriction is defined by the param- 
eters / and w as can be seen in the figure. The external 
magnetic field Bo = Bq e z is oriented perpendicular to 
the two-dimensional structure. 

Since we assume the Fermi surface to be cylindrical 
with the symmetry axis parallel to the z direction, wc 
can restrict our calculations to the xy plane. The results 
are applicable to films of arbitrary thickness with an ap- 
propriate value for the London penetration depth [2l|. 
Our results correspond to the symmetry plane at half 
the extension of the sample in the z direction and, ac- 
cordingly, magnetic fields are oriented parallel to e 2 . 



f(r(x),k F ,ie n ) = 
g(r(x),k F ,ie n ) = -i 



2a{x) 



I + a(x)b(x) 
1 — a(x)b(x) 
1 + a{x)b{x) 



(3) 
(4) 



The Riccati amplitudes a(x) and b(x) are solutions to the 
Riccati differential equations: 

hv F d x a{x) + [2e n + I^{r{x))a{x)]a{x) - A(r(x)) = 

(5a) 

hv F d x b(x) - [2i n + A(r{x)) b(x)} b(x) + A^(f(x)) = 

(5b) 

The Riccati differential equations have to be solved along 
real space trajectories r{x) pointing in the direction of the 
Fermi wave vector v F (k F ) using the modified Matsubara 
frequencies e„: 



ie n — ie n H — v F ■ A{x) 
c 

e n = (2n+l)nk B T 



(6) 
(7) 



The function A(f(x)) = |A(f(a;))| e 1 ^^ in the Riccati 
equations is the complex order parameter. The gradi- 
ent of the order parameter phase V(^(r) together with 
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the magnetic vector potential A{r) consitute the gauge- 
invariant superfluid velocity v s {r): 



v s (r) 



1 

2m 



2e 



A(r) 



(8) 



In the following subsection, we will describe the pro- 
cedure how we calculate the magnetic vector potential 
A(f). The distribution of the order parameter phase (f>(r) 
follows from the solution of the gap equation which, to- 
gether with the current equation (|2|), minimizes free en- 
ergy. Thus, the way we calculate the magnetic vector 
potential makes use of the freedom of gauge, but after 
choosing a certain way to calculate A{f), the order pa- 
rameter phase is uniquely defined by the self-consistent 
solution of the gap and the current equation. 

To obtain stable solutions, the Riccati equation for 
a(x) has to be integrated along vp whereas the Riccati 
equation for b(x) in the opposite direction. Appropriate 
initial values for the integration of the Riccati differential 
equations can be constructed by dropping the gradient 
terms. For e„ > 0, this yields 



o(— oo) 
6(+oo) 



A(-oo) 



e n + ^ + |A(-oo)P 
At(+oo) 

£„ + V^ + I A (+°°)I 2 



(9a) 
(9b) 



For the evaluation of the right hand side of the gap and 
the current equation, the components of the quasiclassi- 
cal propagator f(f{x),k,F,ie n ) and g(f(x),kF,ie n ) con- 
structed with the Riccati amplitudes according to equa- 
tions ([3]) and (H]) have to be averaged over the Fermi 
surface of the superconducting material (denoted with 
(■■■}fs)- As stated above, we assume the Fermi sur- 
face to be cylindrical with the symmetry axis along the 
z direction. Finally, these averaged quantities have to be 
summed over all Matsubara frequencies e n smaller than 
the cutoff u) c . For our calculations, we use a cutoff of 
lu c = 50 ksT c which is sufficiently large for the results to 
be independent of uj c . 

Apart from the left and right end of the channel sec- 
tion, all other boundaries of the geometry are assumed to 
be impenetrable. Thus, we can apply the general bound- 
ary conditions for the Riccati amplitudes given by She- 
lankov and Ozana [25| in a very simple form. Vanishing 
transmissibility leads to specular reflection of the trajec- 
tories (outgoing angle is equal to incident angle) and to 
the conservation of the Riccati amplitudes. 
The optical behaviour of the trajectories makes it nec- 
essary to implement a two-dimensional raytracing pro- 
cedure. Depending on position and angle, this leads to 
multiple reflections of the trajectories within the geome- 
try. A similar raytracing procedure for the Riccati equa- 
tions has been used in [26[ to calculate the local density 
of states at polygonal boundaries of d-wave supercon- 
ductors. In our case, the trajectories can be terminated 
either if the left or right channel end have been reached 



or if the length of the trajectory is sufficient for the re- 
sults of the integration to be independent of the initial 
values. 

To solve the above equations, we discretize the area of 
the channel section and solve the equations at the grid 
nodes. Numerical costs can be reduced significantly by a 
local variation of the grid width. In the direct vicinity of 
the constriction, a rather fine grid is needed whereas less 
nodes are sufficient in the channel. 



B. Magnetic Fields 

To account for magnetic fields, the two self-consistency 
equations (JTJ) and @ have to be completed by the 
Maxwell equation 



V x B 



-in 



■3s 



and the definition of the magnetic vector potential 



V x A = B 



(10) 



(11) 



The total magnetic vector potential has two contribu- 
tions. First, there is the contribution A$ which represents 
the externally applied magnetic field. Second, the super- 
currents flowing in the junction cause a magnetic vector 
potential A c . Combining both yields the total magnetic 
vector potential: 



A = A, 



A n 



(12) 



Choosing Coulomb gauge (VA = 0), it follows from |TT 
and (fTTjl for the contribution evoked by the supercurrents 
that 



AA r 



4-7T 



-Js 



(13) 



In order to solve this Poisson equation for the magnetic 
vector potential A c , we use the two-dimensional Green's 
function of the Laplacian. This leads to 



Ac{r) = - 



4tt 



1 

27 



In |f- 7*1 •£(/)• dV (14) 



Free space boundary conditions for r — > oo are con- 
tained in the Green's function whereas the behaviour 
in the vicinity of the sample is determined by the self- 
consistency conditions |T|) and |J2J), see subsection IIII Dl 
In a two-dimensional discretization of the geometry, this 
integral can easily be solved analytically for each grid cell 
and then summed over all cells. 

In order to have consistent behaviour of the magnetic 
vector potential at the connection to the long channel, 
we assume the currents to be constant away from the 
constriction (for |x| > L) and also integrate over these in 
equation (fl4|) . 

After calculating the magnetic vector potential A c (r) 
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from the supercurrents via eq. (fl4|) , we make use of the 
freedom of gauge. The phase gradients of the order pa- 
rameter V</>(r) together with the total magnetic vector 
potential A(r) lead to the final form of the supercurrents 
(see eq. ©). This implies that the resulting phase gradi- 
ents depend on the form of the magnetic vector potential 
and vice versa. To minimize both the phase gradients 
and the magnetic vector potential, we shift the vector 
potential A c (r) to minimal absolute values. This can 
be done via a simple homogeneous offset which neither 
changes the corresponding magnetic field B — V x A nor 
conflicts with the Poisson equation (fl~3f . Minimizing the 
absolute values of the phase gradient and the magnetic 
vector potential leads to more robust convergence of the 
iterative procedure which we use to solve the gap and 
the current equation. This iterative procedure will be 
described in detail in subsection IIII Dl 

To complete the total magnetic vector potential, we 
include external magnetic fields Bq = Bq e z via Aq = 
—Boy e x which complies with the Coulomb gauge chosen 
for the magnetic vector potential A and which minimizes 
absolute values. 



C. Boundary Conditions for the Order Parameter 

In the problem considered in the present work, two 
factors leading to pair breaking are relevant. First, a 
phase difference is applied which leads to a transport 
current and thus suppresses superconductivity. Second, 
the external magnetic field leads to screening currents 
and pair breaking. In order to apply correct boundary 
conditions at the connection to the channel, these two 
effects have to be clearly separated. 

Since we study Josephson junction behaviour, the max- 
imum transport current will be given by the Josephson 
critical current which is much lower than the critical cur- 
rent of a massive superconducting lead. Thus, the trans- 
port current will not significantly lower the gap ampli- 
tude in the channel, far away from the constriction. For 
vanishing external magnetic field, it is therefore a very 
good approximation to use the temperature dependent 
bulk gap amplitude A 00 (T) as boundary condition at the 
left and right end of the section of the channel enclos- 
ing the constriction. Furthermore, for vanishing external 
magnetic field, the phase is homogeneous over the cross 
section of the channel and a phase difference can easily 
be applied by setting the phase boundary conditions to 
fl , L = ±A0/2. 

The situation changes in the presence of an external 
magnetic field. The magnetic field leads to screening 
currents and lowers the gap amplitude in the channel. 
To gain correct behaviour for higher magnetic fields, it 
is thus essential to include the lowering of the gap am- 
plitude in the channel. With an external magnetic field, 
we thus proceed in two steps. First, we self-consistently 
calculate the amplitude and phase of the order parameter 



for zero transport current (corresponding to zero phase 
difference) without any boundary conditions at the left 
and right section ends. This will give us the correct self- 
consistent solution with screening but without transport 
currents and allows for a lowering of the gap amplitude as 
well as for variations of the phase over the cross section of 
the channel. In the second step, we use the amplitude and 
phase of the order parameter found at the left and right 
section end in the first step as boundary conditions for 
finite transport current. Therefore, the phase has to be 
shifted by —Atfi/2 at the left section end and by +A0/2 
at the right section end whereas the amplitude has to 
be retained. Herewith, the pair-breaking influence of the 
external magnetic field and the screening currents is in- 
cluded for all amplitudes of the transport current. This 
two-step procedure has to be repeated for every value of 
the magnetic field. 

The gradients of the order parameter phase V0(r) to- 
gether with the magnetic vector potential A{r) consti- 
tute the gauge-invariant superfluid velocity v s (r) given 
in eq. ([5]). The current-phase relation thus has to be ex- 
pressed in terms of the gauge-invariant phase difference 
7 with 

7 = A0- — / Adl (15) 

where $o = is the nux quantum and L, R indicate the 
left and right section end. 

D. Iterative Procedure 

Let us now comment on how we gain the self-consistent 
solutions for the order parameter A(r), the supercurrents 
] s {r) and the magnetic vector potential A(r). The self- 
consistency equations ([TJ and ^ serve as iteration rules 
with which we improve an initial guess on a discrete grid. 
In each iterative step, we first have to solve the Riccati 
equations ([5]) for all Matsubara frequencies < e n < uj c 
and all Fermi velocities vf at every grid node. Second, we 
construct f{r{x) 1 kF,ie n ) and g{r{x),kF,i£ n ) from the 
solutions for a{x) and b(x) via ^ and and plug these 
into the gap ([1]) and current equation ([2]). Third, we cal- 
culate the magnetic vector potential A c (f) caused by the 
supercurrents via (|14|) and construct the total magnetic 
vector potential A according to (fT2"|) . With the new con- 
figurations for A(r), f s (f) and A(f), we can attempt the 
next iterative cycle. 

This iterative procedure has to be repeated until con- 
vergence. Two criteria can be utilized to determine the 
level of self-consistency of the solution. It is straightfor- 
ward to observe the change of the functions A(r), j s (f) 
and A(r) in every iterative step and to continue until it 
is sufficiently small. Since the rate of convergence in a 
single iteration can be very small, a more significant ap- 
proach is to verify current conservation which provides a 
stringent test of self-consistency [13, Hl| • 
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In the formulation of the microscopic Eilenberger equa- 
tions used in the present work, the order parameter in the 
gap equation retains a complex quantity. Thus, iterating 
the gap equation yields both the amplitude and the phase 
distribution of the order parameter self-consistently. The 
magnetic vector potential is calculated from the currents 
and the external magnetic field in every iteration and the 
system of equations is thus solved entirely self-consistent 
including the magnetic vector potential. 

At first sight, the free space Green's function (fT4")l for 
the magnetic vector potential does not provide bound- 
ary conditions for the surface of the geometry. Since the 
Green's function (jT4j) is embedded in a self-consistent cal- 
culation, the compliance with the boundary conditions is 
ensured by the gap (p} and the current equation ([2]). A 
self-consistent solution of the Eilenberger equations au- 
tomatically leads to Vj s = and respect of the boundary 
conditions at the surface of the superconductor (n-? s = 
with the surface normal n). Since, in the self-consistent 
calculation, the magnetic vector potential depends on the 
currents and vice versa, the correct behaviour of the mag- 
netic vector potential including VA = follows. Embed- 
ding the Green's function (fT4")) in the self-consistent cal- 
culation thus liberates us from the explicit formulation of 
boundary conditions for the magnetic vector potential. 

It should be noted that an iterative algorithm usually 
produces stable solutions. As discussed in [8J], unstable 
solutions corresponding to saddle points of the free en- 
ergy can only be reached with more elaborate methods. 
However, since we are only interested in stable solutions, 
an iterative approach is controllable and straightforward. 



E. Local Density of States 

Once a self-consistent configuration of the order pa- 
rameter, the supercurrents and the magnetic vector po- 
tential has been found, the normalized local density of 
states follows via 



Applying these normalizations, equation ((T3|) becomes 



N(E, r) 



(Im g(r,k F ,is n -> E + iS) ) „ (16) 



with the definitions of g and (■ ■ - )fs introduced above Q. 



AA C = 



(17) 



Thus, it is obvious that for reasonably large values of 
K, the magnetic vector potential evoked by the currents 
can be neglected and electromagnetic self-consistency be- 
comes unimportant. Nevertheless, for smaller values of n, 
magnetic fields become increasingly important and have 
to be included in the calculations. 



IV. COMPARISON WITH EARLIER RESULTS 

As a first test of our approach and our numerical 
code, we compare our results with those of Ginzburg- 
Landau theory presented in chapter 4.4 of Tinkham's 
textbook [2{|. There, one can find results for the critical 
current of a thin wire or film. In order to compare our re- 
sults with these calculations, we remove the constriction 
(I = 0, w = W/2) and consider temperatures close to T c . 
Additionally, we have to change our boundary conditions 
slightly in order to allow for suppression of superconduc- 
tivity in the channel at high transport current densities. 
This can easily be done via periodic boundary conditions 
at the left and right end of the channel section. After 
these minor changes, our program is ready to consider 
a long thin wire or film and our results perfectly match 
those from Ginzburg-Landau theory. The same results 
for the critical current can also be achieved by a simple 
Doppler shift calculation [30| . i.e. solving the bulk gap 
and current equation with a homogeneous phase gradient 



■ v F ■ V(f>(x) 



(18) 



Normalization 



This consideration provides a second verification. 

As a third test, we compare our results with those for a 
ballistic superconducting microconstriction obtained by 
Zareyan et al. [lfj]. In their work, Zareyan et al. analyt- 
ically solve quasiclassical Eilenberger theory for a ballis- 
tic narrow superconducting channel connecting two bulk 
electrodes for temperatures close to T c . If we remove the 
constriction from our geometry (/ = 0, w = W/2) and use 
our standard boundary conditions for the left and right 
section end described in section [TTTl we can numerically 
reproduce the results given by Zareyan et al. 



All quantities with the dimension of an energy will be 
normalized to ksT c . The characteristic lengthscale we 
use to normalize all lengths is £0 == A Trbo) ■ Total 
currents (integrated over the cross-section of the chan- 
nel) will be given in terms of j tot = gvT^fefiJ^ with 
v F = \v F \. 

With k = \l/£,o and the London penetration depth 
Xf = 4ir(e/c) 2 N(0)v 2 F , we use B = $o/(1 .f 64 ^ } and 
A- ^ 

71 * /(1.7647T 



V. THE JOSEPHSON EFFECT IN SCS 
JUNCTIONS 

In this section, we show results without magnetic fields 
in order to explain in detail how the Josephson effect in 
a ScS junction emerges. As mentioned above, neglecting 
the magnetic fields caused by the supercurrents corre- 
sponds to k — 00. 

In order to explain the Josephson effect in a ScS junc- 
tion in detail, we fix the geometric parameters (I — w — 
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0.05 f ) as well as the temperature (T = 0.5 T c ). In Fig.H 
we show the amplitude (a) and the phase (b) of the order 
parameter for different values of the phase difference as 
well as the current-phase relation (c). Additionally, we 
show the LDOS in the center of the constriction in (d). 

With increasing phase difference, increasing supercur- 
rents flow across the junction. The current concentration 
in the constriction leads to high local current densities 
and causes a suppression of the order parameter ampli- 
tude. This suppression of the order parameter amplitude 
in turn limits the current that can flow across the junc- 
tion. At the same time, the suppression of the order 
parameter amplitude leads to weak coupling between the 
left and right electrode. The suppression of the order pa- 
rameter amplitude caused by the current concentration 
in the constriction thus evokes the Josephson effect in 
ScS junctions. 

As long as the phase difference and the current density 
are small, the suppression of the order parameter ampli- 
tude is negligible and the self-consistent phase distribu- 
tion resembles the solution of the Laplace equation from 
London theory, i.e. from local electrodynamics. The 
stronger the currents, the stronger is the suppression of 
the order parameter amplitude in the constriction and 
the stronger is the deviation of the phase distribution 
from the solution of the Laplace equation. At phase dif- 
ferences close to 7 = 7r, the whole phase drop occurs in 
the constriction and only very small currents flow. 
For the parameters chosen here, we find that the current- 
phase relation strongly deviates from the standard sinu- 
soidal form. As long as the suppression of the order pa- 
rameter in the constriction is small, the current-phase 
relation is basically linear. However, as the suppres- 
sion of the order parameter becomes stronger, the cur- 
rent reaches the Josephson critical current and vanishes 
at 7 = 7r. In Fig. [31 we show the current flow pattern 
in the vicinity of the constriction at the critical current 
(7 = 0.75 tt). 

In Fig. [3 (d), we show the LDOS in the center of 
the constriction. With increasing phase difference, pro- 
nounced Andreev bound states occur at energies below 
the gap (\E\ < A oc (T)). At intermediate values of the 
phase difference, additional structures in the LDOS ap- 
pear close to the main peak of the bound state which 
stem from multiple reflections in the constriction. These 
substructures will be discussed in detail in the next sec- 
tion. At the maximum phase difference of 7 = tt, the 
Andreev bound states occur at an energy of E = and 
exhibit the highest peak value. 

By considering the amplitude and the phase of the 
order parameter, the currents and the local density of 
states, we gain a very complete picture of the Joseph- 
son effect in ScS junctions. We can state that the su- 
percurrent across the constriction-type weak link is in- 
deed mediated by the Andreev bound states. Cooper 
pairs incident on one side of the junction are transformed 
into quasiparticles that travel across the junction phase- 
coherently and recompose to Cooper pairs on the other 
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Figure 2: Josephson effect in a ScS junction with I = w — 
0.05 £0 and k = 00 at a temperature of T — 0.5 T c . In (a) 
and (b), the order parameter amplitude and phase are plot- 
ted along the x axis for several values of the phase difference 
across the junction. The constriction length is indicated with 
the vertical dashed lines. In (c), we show the corresponding 
current-phase relation with the Josephson critical current in- 
dicated by the dashed lines. The LDOS in the center of the 
constriction (x = 0, y = 0) is plotted in (d). 
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Figure 3: Current flow pattern in the vicinity of the constric- 
tion for I — w — 0.05 £o> K — oo and T = 0.5 T c at the critical 
current (7 = 0.75 71"). For better visibility, we plot only every 
fourth point out of the grid used in our calculations. 



side. 

In order to explain the relevance of self-consistency, we 
compare the self-consistent results with those of a S-N-S 
model for the order parameter: 



A(r) 



j(T ). e +iA0/2 



. x < — I 
,\x\<l 

, X > +1 



(19) 



With this model, we calculate the supercurrents for 
phase differences < 7 < n and extract the current- 
phase relation and calculate the LDOS in the center of 
the constriction (see Fig.|4j). The parameters used here 
are the same as for the self-consistent solutions in Fig. [5] 
When using the S-N-S model, one implicitly assumes that 
the junction exhibits weak link behaviour. This justifies 
the suppression of the order parameter and allows for a 
rapid variation of the phase. The result is a more sinu- 
soidal current-phase relation because the full phase drop 
occurs within the constriction (see Fig.Q](a)). In the self- 
consistent solution however, part of the phase variation 
occurs in the channel and the effective phase difference in 
the direct vicinity of the constriction is smaller, at least 
for small total phase differences (cf. Fig. [5] (b)). 
The position of the bound states for the S-N-S model 
reflects the fact that the current-phase relation is more 
sinusoidal (see Fig.Q](b)). Additionally, it is remarkable 
that the bound states in the S-N-S model do not exhibit 
additional substructure. 



VI. INFLUENCE OF TEMPERATURE AND 
GEOMETRIC PARAMETERS 

In the first part of this section, we analyze the tempera- 
ture dependence of the Josephson effect in ScS junctions. 
In the second part, we will focus on the influence of the 
geometric parameters that define the constriction. 

In order to investigate the temperature dependence of 
the Josephson effect in ScS junctions, we keep the geo- 
metric parameters fixed (I = w = 0.05 £0) and neither 
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Figure 4: Current-phase relation (a) and LDOS in the center 
of the constriction (b) for the non-self-consistent S-N-S model 
given in eq. (|19p . The self-consistent current-phase relation 
from Fig. [2] (c) is also plotted in (a) for comparison. For this 
figure, we use the same parameters as for Fig. [5] 



consider external magnetic fields nor those caused by the 
supercurrents {Bq = 0, K = 00). 

To study the influence of the temperature, we calculate 
the order parameter A(r), the currents j[r) and the local 
density of states N(E) for different values of the phase 
difference 7 at temperatures from T = 0.1 T c to 0.9 T c . 
From the resulting current configurations, we extract the 
current-phase relations which we show in Fig. [3] (a). In 
Fig. [5] (b), we plot the amplitude of the order parame- 
ter in the center of the constriction for the same values 
of T and 7 as for the current-phase relations. For the 
two highest temperatures, we increase the cutoff w c used 
in the gap equation and the current equation in order 
to include a sufficient number of Matsubara frequencies 
[w c = 100 k B T c for T = 0.7 T c and u c = 200 k B T c for 
T = 0.9 T c respectively). 

At temperatures close to T c , the current-phase rela- 
tion is highly sinusoidal. The lower the temperature, 
the stronger are the deviations of the current-phase re- 
lation from sin(7/2). For the geometry used here, the 
current-phase relation becomes multivalued at temper- 
atures T < 0.3 T c . As stated above, our iterative ap- 
proach yields only the stable branch of the solution. The 
unstable branch of the current-phase relation is roughly 
sketched in the diagram (dashed lines, see e.g. 0,11]). 
The multivalued character of the current-phase relation 
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Figure 5: Current-phase relations for I — w — 0.05 and 
ft = oo at different temperatures are shown in (a) and the 
corresponding amplitude of the order parameter in the center 
of the constriction in (b). Points with solid guidelines are 
the calculated stable branches whereas the dashed lines are a 
rough sketch of the unstable ones. Inset in (a): Temperature- 
dependence of the critical current (points with solid guideline) 
and the Kulik-Omel'yanchuk (dashed) and the Ambegaokar- 
BaratofF results (dash-dotted) for comparison. 



implies that the current increases even for phase differ- 
ences 7 > 7r and marks the transition from Josephson 
junction behaviour to bulk current flow. 
In the inset of Fig. [5] (a), we plot the dependence of the 
critical current on temperature (points with solid guide- 
line). For comparison, we plot the Kulik-Omel'yanchuk 
result from [31| for a point contact in the clean limit 
(dashed line) as well as the Ambegaokar-Baratoff result 
for a tunnel junction (dash-dotted) from 32]. Close to 
T c , the critical current exhibits a linear decrease with 
increasing temperature and the self-consistent results of 
our calculations coincide with those for the point con- 
tact. At lower temperatures however, the self-consistent 
treatment leads to a stronger increase of the critical 
current than predicted for the point contact. In the 
Kulik-Omel'yanchuk results for the point contact, the 
current vanishes at 7 = 7r for all temperatures whereas, 
in the self-consistent treatment, the currents increase for 
7 > 7t 0, HI- At low temperatures, this leads to higher 
values of the critical current in the self-consistent calcu- 
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Figure 6: Local density of states in the center of the constric- 
tion at three different temperatures (T = 0.1 T c , 0.5 T c and 
0.9 T c from top to bottom) for different values of the phase 
difference 7 as indicated in the diagrams. Inset: Position of 
the Andreev bound states at T — 0.9 T c as a function of the 
phase difference 7 (dots) and E = Aoo(T) cos(7/2) for com- 
parison. For this figure, we use the same parameters as for 
Fig.0 



lation. In calculations for a more point contact-like con- 
striction (I = w = 0.01£o)) we find that this correction 
due to self-consistency remains. 

In Fig. [HI we show the local density of states in the 
center of the constriction for three different tempera- 
tures. With increasing phase difference and thus increas- 
ing transport current, Andreev bound states occur in the 
constriction. The larger the phase difference, the lower 
is the energy of the bound states. At T = 0.1 T c , the 
Andreev bound state for 7 = 7r still has finite energy. At 
higher temperatures, when the current-phase relation is 
not multivalued, the bound states have zero energy for 
the maximum phase difference 7 = it. 
For T — 0.9 T c , we additionally plot the position of the 
Andreev bound states as a function of the phase differ- 
ence 7 (inset in lowermost panel of Fig. [5]). At temper- 
atures close to the critical temperature T c , the position 
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Figure 7: Current-phase relation (a) and order parameter am- 
plitude in the center of the constriction (b) for different values 
of the constriction length I as indicated in the diagrams. For 
this figure, we use T — 0.5 T c , k — oo and w = 0.01 £o and do 
not consider magnetic fields. Lines are guides for the eye. 

of the bound states is given by E = ±A OQ (T) cos(7/2) 
which is plotted in the inset for reference. At lower tem- 
peratures, deviations occur as can be seen in the upper 
two panels of Fig. [HI 

Primarily at low temperatures, the Andreev bound states 
do not only have one peak but exhibit additional local 
maxima at energies above the main peak. This substruc- 
ture is due to multiple reflections of quasiparticles within 
the constriction. If we consider the center of the constric- 
tion and measure the angle of the trajectories a with re- 
spect to the x axis, an additional reflection within the 
constriction sets in at angles a n with 

w 

tan(a„) = (2n+l)-— , n = 0,1,2,... 

The intervals between the angles a n each correspond to a 
local maximum of the bound state. The effective distance 
over which the phase drop occurs along the trajectory in- 
creases with every additional reflection and the effective 
phase gradient thus decreases. 

In the case of the S-N-S model, the substructure of the 
bound states does not occur because the phase drop oc- 
curs in the direct vicinity of the constriction. Reflections 



Figure 8: Current-phase relation (a) and order parameter am- 
plitude in the center of the constriction (b) for different values 
of the constriction width w as indicated. The dashed lines cor- 
respond to a channel without constriction. For this figure, we 
use T = 0.5 T c , k = oo and / = 0.01 and do not consider 
magnetic fields. Lines are guides for the eye. 



of the trajectories do not change the effective distance 
over which the phase drop occurs along the trajectories. 

Now let us address the influence of the geometric pa- 
rameters on the Joscphson effect. 

In Fig. we plot the current-phase relation and the order 
parameter amplitude in the center of the constriction for 
different values of the constriction length I. The temper- 
ature is set to T — 0.5 T c , the width of the constriction 
to w — 0.01 £o &ud we do not consider magnetic fields. 
If the constriction is shorter, more current flows and the 
suppression of the order parameter amplitude is stronger. 
Thus, for shorter constrictions, the current-phase relation 
is more sinusoidal and the Josephson critical current is 
higher than for longer ones. The decrease of the critical 
current with increasing length of the constriction is con- 
sistent with the results found by Zareyan et al. [l(| for 
temperatures close to T c . 

In Fig. [51 we show the current-phase relation and the 
order parameter amplitude in the center of the constric- 
tion for different values of the constriction width w. As 
before, the temperature is set to T — 0.5 T c , the length 
of the constriction to I = 0.01 £o an d we do not consider 
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magnetic fields. For reference, we show the current-phase 
relation without a constriction but with periodic bound- 
ary conditions as described in section IIVI (dashed lines 
in both Fig. [5] (a) and (b)). The dashed lines thus corre- 
spond to a long superconducting lead with homogeneous 
current flow or equally a homogenous bulk current. 

As expected, the width of the constriction strongly in- 
fluences the Josephson critical current. The narrower the 
constriction, the stronger is the current concentration in 
the constriction and the stronger is the suppression of the 
order parameter amplitude. Hence, less current flows at 
the same phase difference for a narrower constriction. 



VII. EFFECTS OF SCREENING 

From this section on, we will include the effects of 
magnetic fields in our considerations. In this section, we 
will focus on the influence of the parameter k = \l/£,o 
whereas in the next section, we will examine the effect of 
an external magnetic field. 

So far, we neither considered external magnetic fields 
nor those caused by the supercurrents. Neglecting the 
latter corresponds to large values of n which can easily 
be seen from equation (fTT)) . In our case, the lateral ex- 
tension of the geometry is W — 1 £o and thus not much 
happens for k > 1. Nevertheless, for smaller values of k, 
the distribution of the currents and the magnetic fields 
changes drastically. 

In Fig. [21 we show the current density at the Josephson 
critical current for k = 1 (a) and k = 0.1 (b). To calculate 
these configurations, we use T = 0.5 T c and I — w — 
0.05 £o- For k = 1, the current density is homogeneous 
over the cross section of the channel. Nevertheless, in 
the vicinity of the constriction, the current is strongly 
concentrated. For k = 0.1, the current is concentrated 
on the surface of the sample. In the channel, the current 
distribution thus changes drastically. However, since the 
constriction size is equal to the London penetration depth 
Xl = «£o with k = 0.1 used here, the current density is 
homogeneous in the constriction. 

In Fig. 1101 we plot the magnetic fields corresponding 
to the current configurations in Fig. [51 The strong con- 
centration of the current on the surface of the sample 
leads to the disappearance of the magnetic field in the 
interior of the channel. Nevertheless, the magnetic field 
in the vicinity of the constriction does not depend much 
on the value of k as long as the dimensions of the con- 
striction are smaller than the London penetration depth 
(Al > l,w). Accordingly, the current-phase relation and 
the dependence of the order parameter amplitude in the 
center of the constriction on the phase difference remain 
unchanged for k > 0.1. Thus, we can conclude that the 
Josephson effect is independent of k as long as Xl > I, w. 

In Fig. [Til we show the local density of states in the 
center of the constriction for T — 0.5 T c and I = w = 
0.05 £o- The inhomogeneous configuration for small val- 
ues of k leads to a splitting of the bound states. For the 



0) 0.5 



ig 



-0.5 

-1 -0.5 0.5 1 




-1 -0.5 0.5 1 

x/fo 



Figure 9: Current density at the critical current for n — 1 (a) 
and k = 0.1 (b). Dark shading represents low current density 
wheres bright areas exhibit high current density. Arrows in- 
dicate the direction of the current flow. For small values of 
K, the current is strongly concentrated on the surface of the 
sample. 



case of k = 0.1 and 7 = 0.77 n, we discuss the origin of 
the peaks labeled 1,2,3 in the upper panel of Fig. [Til in 
detail. 

Peak 1 originates from trajectories which are oriented 
parallel to the x axis and thus pass through regions with 
small current density. Peaks 2 and 3 stem from trajecto- 
ries which are reflected into the lower (y < 0) and upper 
(y > 0) part of the channel, i.e. in regions with high 
current density. Thus, the concentration of the super- 
currents on the surface of the sample leads to differing 
contributions of the trajectories for different angles and 
yields a splitting of the bound state. 



VIII. INFLUENCE OF THE EXTERNAL 
MAGNETIC FIELD 

In this section, we will finally examine the influence 
of an external magnetic field on the Josephson effect in 
a ScS junction. Therefore, we set the temperature to 
T = 0.5 T c , the geometric parameters to I = w — 0.05 £0, 
the parameter k = 1 and apply an external magnetic field 

Bo = Boi z . At the end of this section, we will discuss the 
temperature dependence of the influence of the external 
magnetic field. 

Before turning to the results of our calculations, let us 
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Figure 10: Magnetic field distribution B — B e z for « = 1 (a) 
and k — 0.1 (b). We multiply the magnetic field B with k 2 in 
order to have equal scaling of the vertical axis. The geometry 
of the sample is indicated with the black lines. For smaller 
values of k, the current is concentrated on the surface of the 
sample and the magnetic field is suppressed in the interior. 

briefly discuss the normalization we use for the magnetic 
field in terms of more familiar quantities. Following the 
procedure presented in [33j , we calculate the bulk upper 
critical field B c2 using Eilenberger theory (34j]. For s- 
wave superconductors with cylindrical Fermi surface and 
with the magnetic field oriented parallel to the cylinder 
axis, we find 

\T C ) J sinhw V J 

At T = 0.5 T c , the normalized upper critical field thus is 
B c2 = 4.51, or B c2 = 0.424- (- -T c ) respectively 

An external magnetic field leads to screening currents 
which flow on the surface of the superconductor. Thus, in 
our case, the currents are a mixture of transport currents 
across the Josephson junction and screening currents. In 
Fig.HH we show the current density in the section en- 
closing the constriction for an external magnetic field of 
Bq = 2 without transport currents (7 = 0, (a)) and with 
critical transport currents (7 = 0.747T, (b)). In Fig. [TBI 
we plot the corresponding magnetic field distributions 
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Figure 11: Local density of states in the center of the constric- 
tion for two different values of n (k — 0.1 in the upper panel, 
k = 1 in the lower panel). For this figure, we use T = 0.5 T c 
and I = w = 0.05 £o- Please note the different scaling of the 
vertical axis for the two panels. 



Be z . 

Without transport currents, we find screening currents 
which flow on the surface of the channel and which lower 
the magnetic field inside the sample. Since the width 
of the channel is equal to the London penetration depth 
(A j, = k£o with k = 1 used here), the magnetic field is 
not completely suppressed. In the constriction, a small 
circulating current flows which is directed opposite to the 
screening currents and which leads to a local enhance- 
ment or antiscreening of the magnetic field. 
With a finite phase difference applied between the ends 
of the section, a transport current flows. This leads to a 
breaking of the fourfold symmetry visible in Fig.[T2] (a) 
and to a combined current configuration consisting of 
screening as well as transport currents. Consequently, 
the magnetic field distribution shows an asymmetry. The 
total magnetic field is decreased for y < and increased 
for y > 0. Far away from the sample, the external mag- 
netic field is recovered, but this happens on a much larger 
scale and can not be seen in Fig. [13] 

In Fig. 1141 we show the current-phase relation and 
the order parameter amplitude in the center of the con- 
striction for different values of the external magnetic 
field. For small values of the external magnetic field, 
the strength of the currents is gradually reduced as is 
the order parameter amplitude in the constriction, but 
the functional form of the current-phase relation and 
the order parameter amplitude is basically unchanged. 
At larger values of the external magnetic field however, 
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Figure 12: Current density with applied external magnetic 
field (k = 1, Bq = 2) without transport current (7 = 0, 
(a)) and with critical transport current (7 = 0.74 n, (b)). 
Again, dark shading signifies low current density whereas 
bright shading signifies high current density. Arrows indicate 
the direction of the current flow. Without transport currents, 
a small circulating current appears in the constriction which 
leads to local antiscreening (white circular arrow in the con- 
striction in (a)). With an applied phase difference, the total 
current distribution becomes asymmetric and consists of both 
screening and transport currents. 



the transport currents become very small and the order 
parameter amplitude in the constriction remains unaf- 
fected. At the same time, the Josephson effect disap- 
pears. 

The reason for the disappearance of the Josephson ef- 
fect is that with increasing external magnetic field, su- 
perconductivity is suppressed in the channel. Since the 
constriction is much smaller than the London penetration 
depth, the screening currents do not flow across the con- 
striction (see Fig. [12] (a)). This leads to pair breaking in 
the channel and lowers the amplitude of the order param- 
eter. The screening currents mainly flow in the channel 
and superconductivity first breaks down here. However, 
if superconductivity is already suppressed in the chan- 
nel, no current can flow across the constriction and the 
Josephson effect does not occur. 

The external magnetic field which leads to suppression 
of the Josephson effect is about twice the approximate 
upper critical field B C 2 derived above. To verify the va- 
lidity of this unusually large value of the magnetic field 
leading to the suppression of superconductivity in the 
channel, we carried out calculations for a longer chan- 
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Figure 13: Magnetic field distribution B = B e 2 for applied 
external magnetic field (k = 1, Bo = 2) without transport 
current (7 = 0, (a)) and with critical transport current (7 = 
0.74 7r, (b)). The geometry of the sample is indicated with 
the black lines. Screening currents lead to a reduction of 
the magnetic field in the interior of the sample. Additional 
transport currents increase the magnetic field on one side of 
the sample and decrease it on the other side. 



nel with the same width (L = 10£o, W = l£o), but 
without a constriction. The amplitude of the order pa- 
rameter in the center of the constriction as well as in the 
center of the long channel for 7 = is plotted in the 
inset in Fig. [14] (b). In the channel, superconductivity is 
completely suppressed at somewhat lower values of the 
external magnetic field, namely at about Bq * 8.5. 

Since our two-dimensional geometry implies transla- 
tional invariance in the z direction, the high critical field 
found in our calculations can be understood in terms of 
the parallel critical field H c ^ [35J, |36[ . Because of the in- 
finite extension of the geometry in the z direction, the 
external magnetic field is oriented parallel to the sur- 
faces. Up to an external magnetic field iJ C 3 = 1.695 H C 2 
oriented parallel to the surface of a superconductor, there 
is a superconducting surface sheath of thickness ~ £ (T) 
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Figure 14: Current-phase relation (a) and order parameter 
amplitude in the center of the constriction (b) for different 
values of the external magnetic field Bo. Inset in (a): De- 
pendence of the Josephson critical current on Bo- Inset in 
(b): Amplitude of the order parameter in the center of the 
constriction (solid line) and in the center of a long channel 
without constriction (dashed line) for 7 = as a function 
of the external magnetic field. For the whole figure, we use 
T — 0.5 To, I = W = 0.05 £0 and k = 1. Lines are guides for 
the eye. 



while A — > in the interior. If the geometry consists of 
two coplanar surfaces with spacing d < £, nucleation is 
possible at even higher values of the external magnetic 
field, growing with ~ 1/d (see [13, HI] and references 
therein). Since the channel width is W = l£o, the value 
of the external magnetic field at which superconductiv- 
ity is suppressed is in good agreement with H c ^. In the 
center of the constriction, the spacing between opposite 
surfaces is even smaller and nucleation occurs at even 
higher values of the external magnetic field. 

From the inset in Fig. [14] (b), we conclude that the 
current across the constriction should be suppressed at 
about Bo ea 8.5 since this value of the external mag- 
netic field leads to suppression of superconductivity in 
the channel. The section used for the calculations with 
a constriction (L = 2£o) is too short to reach the undis- 
turbed channel, far enough from the weak link. However, 
using a longer section would lead to a badly defined phase 
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Figure 15: Local density of states in the center of the con- 
striction for three different values of the external magnetic 
field (Bo = 0, 4 and 8 from bottom to top). For this figure, 
parameters are set corresponding to Fig. 1141 Please note the 
different scaling of the vertical axis for the three panels. 
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Figure 16: Temperature dependence of the critical current 
without magnetic fields (k — 00, Bo — 0, solid guideline) 
and with external magnetic field (k — 1, Bo as indicated) 
for I = w = 0.05 £o- For comparison, we plot the Kulik- 
Omel'yanchuk result (|3lj. dashed line) and the Ambegaokar- 
Baratoff result ( 32], dash-dotted). 



14 



difference across the junction and to excessive numerical 
costs. Our main conclusion, the fact that the external 
magnetic field primarily affects the channel and not the 
constriction, remains unchanged, however. 
In a thin film with small extension in the z direction, the 
relevant critical field for the suppression of superconduc- 
tivity is the upper critical field H C 2- The occurrence of 
the parallel critical field H c s in our results stems from the 
assumption of translational invariance in the z direction. 

In Fig.[TSl we plot the local density of states in the 
center of the constriction for three different values of the 
magnetic field. At intermediate strength of the exter- 
nal magnetic field, we find a splitting of the Andreev 
bound states, but the main features of the LDOS with- 
out a magnetic field persist. We still find the energy gap 
and pronounced bound states if a phase difference 7 ^ 
is applied. For 7 = 7T, the bound states have zero en- 
ergy which indicates Josephson junction behaviour. For 
very strong external magnetic fields, the local density of 
states is strongly disturbed. Even without a phase differ- 
ence, the energy gap is partly filled. A phase difference 
still leads to bound states, but their spectral weight is 
strongly reduced. For 7 = 7r, the bound states have fi- 
nite energy and the Josephson effect is absent. 

In Fig.[16l we show the temperature dependence of the 
critical current for I = w = 0.05 £0 without (k = 00, 
Bq = 0) and with external magnetic field (k = 1, Bq 
as indicated in the figure). For comparison, we plot 
the Kulik-Omel'yanchuk result for a point contact in the 
clean limit [3l| and the Amgebaokar-Baratoff result for 
a tunnel junction [32|. We find that the universality of 
the inclination of j c {T) close to T c persists even in an 
external magnetic field. 

IX. CONCLUSIONS 

We present self-consistent solutions of microscopic 
Eilenberger theory for a two-dimensional model of a ScS 
Josephson junction. Magnetic fields, external ones as well 



as those evoked by the supercurrents, have been included 
and the relevant equations have been solved numerically 
without further assumptions. 

In the self-consistent calculation, the Josephson effect 
appears without further input. Unlike in a non-self- 
consistent calculation, we do not have to make assump- 
tions about the existence and the behaviour of the weak 
link. By considering the self-consistent results, we ex- 
plain in detail how the Josephson effect in a ScS junc- 
tion emerges. Taking into account the order parameter 
amplitude and phase, the currents and the local den- 
sity of states, we point out that the Josephson effect is 
a result of phase-coherent quasiparticle transport across 
the junction. Localized quasiparticle states (Andreev 
bound states) appear in the junction if a supercurrent 
flows. These bound states, together with the correspond- 
ing current-phase relations, have been calculated and dis- 
cussed in detail for different geometric parameters, dif- 
ferent temperatures, different values of k = and 
different values of the external magnetic field. 

We show that the Josephson effect in ScS junctions 
is independent of k as long as the London penetration 
depth is larger than or comparable to the size of the con- 
striction. We also show that the Josephson effect in ScS 
junctions is very insensitive to external magnetic fields. 
External magnetic fields mainly influence the channel 
leading to the weak link and suppress superconductiv- 
ity there. The Josephson effect persists up to the upper 
critical field of the channel. 

In contrast to the non-self-consistent results, the self- 
consistently calculated local density of states exhibits An- 
dreev bound states with peculiar substructure due to the 
geometry and magnetic fields. The local density of states 
should be experimentally accessible via STM/STS. 
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